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Abstract. We present the first measurement of the large-scale cross-correlation of Lya 
forest absorption and Damped Lyman a systems (DLA), using the 9th Data Release of the 
Baryon Oscillation Spectroscopic Survey (BOSS). The cross-correlation is clearly detected on 
scales up to 40 h~ l Mpc and is well fitted by the linear theory prediction of the standard Cold 
Dark Matter model of structure formation with the expected redshift distortions, confirming 
its origin in the gravitational evolution of structure. The amplitude of the DLA-Lya cross- 
correlation depends on only one free parameter, the bias factor of the DLA systems, once the 
Lya forest bias factors are known from independent Lya forest correlation measurements. 
We measure the DLA bias factor to be bp = (2.17±0.20)/3^ 22 , where the Lya forest redshift 
distortion parameter /3p is expected to be above unity. This bias factor implies a typical 
host halo mass for DLAs that is much larger than expected in present DLA models, and is 
reproduced if the DLA cross section scales with halo mass as Mfi, with a = 1.1 ± 0.1 for 
f3p = 1. Matching the observed DLA bias factor and rate of incidence requires that atomic 
gas remains extended in massive halos over larger areas than predicted in present simulations 
of galaxy formation, with typical DLA proper sizes larger than 20 kpc in host halos of masses 
~ 10 12 M & . We infer that typical galaxies at z ~ 2 to 3 are surrounded by systems of atomic 
clouds that are much more extended than the luminous parts of galaxies and contain ~ 10% 
of the baryons in the host halo. 

Keywords: cosmology: large-scale structure — cosmology: intergalactic medium — cosmol- 
ogy: galaxy formation — quasars: absorption systems 
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1 Introduction 

Damped Lya systems, denned as absorption systems with neutral hydrogen column density 
Nhi > 2 x 10 20 cm -2 (hereafter referred to as DLAs), are a powerful probe of the physical 
evolution of gas that has condensed to high density to become self-shielded and atomic, and is 
presumably in the process of forming galaxies. Spectroscopic surveys of quasars to search for 
DLAs are usually performed at z > 2, because the absorption spectra can be observed from 
the ground at wavelengths longer than the atmospheric cutoff. These surveys have revealed 
the rate of incidence of DLAs to be ~ 0.2 per unit of redshift at z = 3, slowly increasing with 
redshift, and their total gas content to be close to Qdla — 10~ 3 , or a few percent of all the 
baryons in the universe (see the review [1]). This baryonic mass of DLAs is about one third 
of all the mass that is contained today in stars [2, 3], and is comparable to the mass of stars 
that had formed in the universe by z = 2 (although the total stellar mass at this redshift is 
substantially uncertain; see [4], for a review). 

These observational facts suggest that DLAs are connected with the gas clouds that 
are responsible for forming galaxies at high redshift. Moreover, they imply that once the 
gas that has gravitationally collapsed into dark matter halos becomes self-shielded against 
ionizing radiation and mostly atomic, it must either remain in atomic form in the DLAs 
for a time comparable to the Hubble time before forming stars, or else be ionized, expelled, 
and continuously replaced by other gas accreting onto halos and recombining. The reason is 
that if atomic gas were to quickly form molecular clouds and stars soon after recombining, 
then its baryonic content would always be much smaller than that of stars during the entire 
epoch of galaxy formation. At the same time, the total cross section for intersecting DLAs is 
much larger than the fraction of the sky covered by the starlight-emitting regions of observed 
galaxies. Therefore, large reservoirs of atomic gas must remain orbiting in halos for long 
periods of time, either covering wide areas around the galaxies where star formation is most 
active, or remaining stable in numerous low-mass halos with little star formation. 

Despite these powerful observational constraints, the precise nature of DLAs remains 
poorly known, primarily because the mass of their host halos and the type of galaxies they 
are associated with has not been observationally determined. The observed rate of incidence 
tells us the product of the number density of DLA systems times their cross section, but 
these two quantities are not separately known. In reality, there can be a diverse population 
of objects with a wide range of cross sections giving rise to the observed DLAs. In the 
context of the Cold Dark Matter model of structure formation by hierarchical merging of 
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halos, we can define a mean proper cross section E(M, z) of a halo of mass M at redshift z 
for producing a DLA in the spectrum of a background source. The observed rate of incidence 
of DLA absorbers per unit redshift, R(z), is given by 



R(z) dz 



c(l + z) 2 
H(z) 



oo 



dMn{M,z)T,(M, z) dz 



(1.1) 



where n(M, z) dM is the comoving number density of halos of mass M within the mass range 
dM, and H(z) is the Hubble constant at redshift z. The observational determination of any 
diagnostic of the halo mass associated with specific DLAs has generally proved extremely 
difficult. One of the possible avenues is to detect the galaxy associated with a DLA absorber. 
Recent progress in the selection strategies and observational techniques has allowed for several 
detections (see [5] and references therein), but the sample size remains very limited owing 
to the faintness of the associated galaxies and the additional difficulty involved in detecting 
them at a very small impact parameter from a bright quasar. 

Another method of characterizing the population of halos that are hosting the DLA 
absorbers is through the large-scale clustering amplitude. In the limit of large scales, any 
population of objects that traces the primordial mass perturbations has a correlation function 
that is proportional to the mass autocorrelation in the linear regime, £ m (r). Hence, the 
autocorrelation of DLA absorbers on large scales in real space is ^o( r ) = b D ^ m (r), where bo 
is the bias factor of DLA absorbers. The bias factor of halos of mass M, bh(M,z), can be 
computed by means of approximate analytic models (e.g., [6]) and accurately predicted with 
numerical simulations [7]. The bias factor of the DLA absorbers is related to that of halos 



In general, the bias factor of halos increases with their mass: halos collapsing out of rare, 
high peaks, with a mass much higher than that of a typical halo, are highly clustered [8]. 
Any measurement of the correlation amplitude of DLAs can be a powerful probe of the 
characteristic halo mass hosting the DLAs: the higher the value of the bias, the more massive 
their typical host halos need to be. 

Measuring the bias factor of DLA absorbers requires either measuring their auto- 
correlation, or their cross-correlation with another tracer population. The principal obstacle 
for measuring the clustering amplitude of DLAs has been the sparseness of quasars at z > 2 
that are bright enough to allow for spectroscopy to detect DLAs in absorption; in addition, 
only ~ 10% of observed quasars yield a detected DLA. The first measurement of the clus- 
tering of DLAs was performed in [9] using the cross-correlation with luminous Lyman break 
galaxies, by measuring redshifts of Lyman break galaxies identified in deep imaging of fields 
around 9 quasars with 11 known DLAs in their absorption spectra. Their result was that 
the bias factor of DLAs is in the range 1.3 < bo < 4. The large uncertainty is due to the 
small size of their DLA sample, which is difficult to increase because of the large observing 
time required to map the area around each DLA. 

The Baryon Oscillations Spectroscopic Survey (BOSS, [10]) in the SDSS-III Collabora- 
tion [11] provides a new opportunity for an accurate measurement of the clustering of DLAs. 
With a tenfold increase of the number of known quasars at z > 2 compared to previous 
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surveys, and a similar increase in the number of detected DLAs, we can attempt to cross- 
correlate the DLAs with any other objects that are measured in the available quasar spectra 
over the same redshift range: these are the DLAs themselves, quasars, the Lya forest, and 
metal line absorbers. Because of the low density of quasars, we can measure these correla- 
tions more easily on large scales, corresponding to the typical separation between neighboring 
BOSS quasars on the sky of 15', or about 15 /i" 1 Mpc (comoving), instead of the smaller 
scales that are probed by deep imaging of Lyman break galaxies in the area adjacent to 
targeted DLAs. Measurements on large scales have the advantage that accurate predictions 
for the correlations can be obtained from linear theory. These large-scale correlations can 
also be measured for metal-line absorbers, and a first detection has already been presented 
by [12]. 

In this paper, we measure the cross-correlation of DLAs with the Lya forest transmit- 
ted flux fraction. This turns out to be the cross-correlation that yields the most accurate 
measurement of the bias factor of DLAs, because of the large number of independent Lya 
forest fluctuations that are probed on every line of sight to a quasar. This cross-correlation 
is proportional to the product of a bias factor for DLAs and a bias factor for the Lya for- 
est. The Lya forest bias factor is independently derived from the observed autocorrelation 
of the Lya forest transmitted flux [13], and therefore the DLA bias factor can be robustly 
inferred from our measurement. We describe the sample of DLAs and Lya forest spectra 
used for our analysis in section 2. The detailed method for measuring the DLA-Lya forest 
cross-correlation is presented in section 3, and the results for the cross-correlation and the 
inferred bias factor of the DLAs are presented in section 4. Finally the results are discussed 
in section 5. 

Throughout this paper we use the flat ACDM cosmology, with Q m = 0.281, = 0.0462, 
h = 0.71, n s = 0.963 and as = 0.8, similar with the best-fit parameters obtained from the 
WMAP analysis in [14]. We note that the values of the bias factors we quote for the Lya 
forest and for DLAs vary in inverse proportion to the assumed value of as- 

2 Data 

This study uses the quasar catalogue of the Data Release 9 (DR9,[15]) of the Baryon Os- 
cillation Spectroscopic Survey (BOSS, [10]), which is part of SDSS-III Collaboration ( [11], 
[16], [17], [18], [19], [20] ). The catalogue, described in detail in [21], contains a total of 
61931 quasars at z > 2.15. The target selection procedure used for identifying the quasar 
candidates for BOSS spectroscopy was presented in [22], and uses the methods described in 
[23], [24], and [25]. 

We now describe the sample of DLAs found in the spectra of these quasars and the set 
of Lya forest spectra that we use for measuring the DLA-Lya cross-correlation. 

2.1 The DLA catalogue 

Our analysis is based on a subset of the DLA catalogue of [26]. Here we briefly summarize 
the DLA detection method and refer the reader to [27] for further details. DLAs are searched 
in each line of sight over a redshift range between z m i n , defined as the redshift where the 
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spectral signal-to-noise ratio per pixel, averaged over a 2000 km s _ window, reaches 2, and 
the quasar redshift. For the purpose of detecting DLAs, the quasar continuum is modeled 
by fitting a modified power-law with smoothly changing index plus Moffat profiles on top 
of the emission lines. The data is correlated with synthetic profiles of increasing column 
densities in order to detect DLA candidates and obtain a first guess of the Nhi- Whenever 
the system has associated metal lines, these are used to further improve the accuracy of the 
DLA redshift. Finally, the column density is obtained by fitting a Voigt-profile to the data. 

The overall DLA catalogue contains a total of 12081 DLA candidates with log Nhi > 20. 
For the purpose of measuring the DLA incidence rate and the total amount of neutral gas 
in the Universe, a statistical sample optimized to achieve high completeness was defined by 
[26], where systematic effects can be quantified using mock data. Here, a high purity of the 
catalogue is our greatest concern to measure the cross-correlation because any inclusion of 
objects that are not real DLAs may systematically decrease the measured amplitude of the 
cross-correlation. Completeness, on the other hand, is less important because eliminating 
a fraction of the real DLAs will increase the error of, but not systematically modify, the 
cross-correlation. 

The first cut we apply is to eliminate DLAs outside the redshift range 2.0 < z < 3.5. 
The few DLAs that are outside this redshift range have very few nearby lines of sight with 
good signal-to-noise ratio in which the cross-correlation can be measured, so we eliminate 
them in order to have a well defined redshift interval of our systems. The standard definition 
of a DLA is an absorption system with a column density Nhi > 10 20 ' 3 cm -2 [28], so strictly 
speaking, the systems we use below this column density are sub-DLAs. We decide, however, 
to adopt a threshold of a factor of two below this definition because (i) systems down to 
10 20 cm -2 are robustly identified and are not expected to sharply change their nature with 
column density and (ii) our measurement of the cross-correlation increases in accuracy with 
the number of systems available. In this paper we refer to all of the systems used for our 
cross-correlation measurement, with Nhi > 10 20 cm" 2 , as DLAs (however, when discussing 
rates of incidence and baryon contents of DLAs later in section 5, these quantities will refer 
to systems with Nhi > 10 20 ' 3 cm -2 ). Systems of even lower column density are less reliably 
detected. We shall see in section 4 that we do not detect any dependence of the clustering 
properties with column density. The number of DLAs that satisfy these criteria is 10512. 

The second cut requires that the continuum-to-noise ratio (CNR) of the quasar spectrum 
is larger than three. The continuum-to-noise ratio is defined as the median value of the ratio 
of the fitted continuum in the DLA detection analysis in [26] to the noise per pixel, over 
the observed Lya forest region for each quasar spectrum. This provides a good estimate of 
the data quality over the region of interest, while being independent of the presence of an 
absorber. We find this criterion to be a good threshold for ensuring the purity of our sample 
without excessively reducing the number of systems. The number of DLAs that survive this 
second cut is 9288. 

The rest of the cuts we apply involve eliminating Broad Absorption Line systems (here- 
after, BALs) that produce broad spectral troughs which, at the low signal-to-noise ratio of 
most of the BOSS spectra, are not easily distinguished from the Voigt profiles of DLAs with 
the superposed Lya forest. Our third cut eliminates all DLAs found in quasars classified as 
BALs in the visual inspection of the DR9 quasars described in [21], leaving 8469 DLAs. In 
addition, in order to remove any BAL contaminants near the Lya emission line that may be 
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Figure 1. Quasar rest- frame wavelength distribution of the DLAs that satisfy the criteria of rcdshift 
and column density range, continuum to noise ratio and velocity separation from the background 
quasar. The red histogram (8189 systems) contains also quasars with BAL systems with a Balnicity 
index Bl < 1000 km s~ , while the blue histogram (7458 systems) contains no systems flagged as 
BAL. Vertical lines show the cut applied to our final sample, 1005 A < A r < 1037 A, reducing to 6780 
the final number of DLAs used in this study. 



too weak to have been identified, our fourth cut eliminates all the DLAs that are within a 
velocity separation v < 5000 km s -1 from the quasar redshift, where 

V _ Z q ~ Z D _ X a - A r 
C 1 + Z q X a 

z q and zd are the redshifts of the quasar and the DLA, \ a is the Lya wavelength, and A r is 
the quasar rest-frame wavelength at which the DLA absorption line is centered. Equivalcntly, 
this condition is A r < 1196 A. Application of this constraint reduces the number of DLAs 
to 7458. We emphasize that most of the 1011 systems eliminated in this velocity range 
are probably real DLAs, but we prefer to eliminate them because they probably contain a 
substantial fraction of BAL absorbers among them. 

Figure 1 shows the distribution of A r for the 7458 DLAs passing all of the above cuts 
as the blue (solid) histogram. The red (dotted) histogram shows the same distribution when 
BALs with a Balnicity Index Bl < 1000 kms" 1 ([29]) are included. The histograms clearly 
show an excess of DLAs in the interval indicated by vertical black lines in figure 1. Moreover, 
quasars with weak BALs are much more common in this interval than at other values of A r . 
This strongly suggests that this excess is due to BAL contamination from the Ly/3 and OVI 
absorption lines, and that additional BALs that may be stronger in OVI than in CIV and 
other detectable lines are probably lurking among the DLAs identified in this interval of A r . 
We therefore apply a fifth cut, eliminating all the DLAs in this interval, with a maximum 
wavelength A r < 1037 A chosen equal to the longest wavelength member of the OVI doublet, 
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Figure 2. Left panel: distribution of the 6780 DLA redshifts. Right panel: neutral hydrogen column 
density distribution. The vertical line on the right panel corresponds to the standard lower limit for 
DLAs, log N m > 20.3. 

and a minimum wavelength A r > 1005 A, corresponding to a velocity v < 6000 km s" 1 to the 
blue of the Ly/3 line. This restriction reduces our DLA sample to 6780 systems. In section 
4, we shall show that the bias factor measured for the DLAs found in this interval is in fact 
much smaller than the ones outside, confirming a large degree of contamination by BAL 
systems. 

The column density and redshift distributions of the final set of DLAs selected for our 
study are plotted in figure 2. The redshift distribution of the DLAs used for our measurement 
of the cross-correlation with the Lya forest peaks at z — 2.2 and declines smoothly as a 
function of redshift. 

2.2 Lya forest spectra 

We select quasars in the redshift range 2.15 < z < 3.5 for the Lya forest spectra to be 
correlated with the DLA positions. The total number of DR9 quasars in this redshift range 
is 58722. We eliminate all the spectra that are identified as BAL, and we also make a 
signal-to-noise ratio cut requiring a median S/N > 0.5 per pixel in the rest-frame wavelength 
interval 1220 A < A r < 1600 A. The number of spectra that are left after these cuts is 52449. 

The co-added spectra released in DR9 are used, removing any pixels that do not pass 
the bit mask and the sky mask as defined in [30]. For this study we define the Lya forest 
region as the pixels with a rest-frame wavelength in the range 1041 A <X r < 1185 A. 

A total of 3047 of these lines of sight contain at least one DLA in the Lya forest region. 
For most of our results we include these lines of sight for our estimate of the cross-correlation 
of DLAs with the Lya absorption, using the correction explained in [30] for removing the 
DLA from the Lya absorption, where the central region of the DLA is masked and the wings 
are corrected using a Voigt profile with the best estimate of the column density. As shown 
in section 4.6, rejecting these lines of sight does not significantly change the results. Since 
DLAs are not easily identified in low S/N spectra, we expect some residual contamination of 
DLAs in our Lya sample. 
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3 Method 



This section describes first the determination of a continuum model for each quasar to infer 
the Lya absorption field. The method we use to compute the cross-correlation of DLAs with 
the Lya forest and its covariance matrix, and to fit a parameterized linear theory model 
to the result, is then presented. We use the co-added quasar spectra of DR9 selected as 
explained in section 2, which contain combined exposures from a single plate after applying 
interpolated sky subtraction and flux calibration corrections [10, 16]. The observed flux fi 
at each pixel i at wavelength Aj = A a (1 + Zi) (where X a is the Lya wavelength and Zi is 
the Lya absorption redshift) is the product of the quasar continuum d and the transmitted 
flux fraction Fi, plus the noise iVj that we approximate as Gaussian with variance (^Vf), as 
provided by the pipeline [16] , 

ft = d F t + Ni = d F( Zi ) [1 + Spi] + Nt , (3.1) 

where F(zi) is the mean value of F as a function of redshift, and 5fi is the Lya transmission 
perturbation. Our analysis of the correlation will use the variable 8fu where often the pixel 
subindex i will be dropped for brevity. 

3.1 Continuum fitting 

We use the Principal Component Analysis (PCA) technique described in [31] for fitting a 
continuum to each quasar. A set of PCA quasar templates with eight components is used to 
do a least-squares fit on each quasar spectrum in the 1216 A < A r < 1600 A region to obtain 
an estimate for the continuum, using equation (3) of [31]. The total number of parameters 
used is 11 (eight for PCA eigenvalues, one for the flux normalization, and two for a redshift 
and mean spectral slope corrections; see Table 1 of [31]). The PCA templates that are used 
span the wavelength interval 1020 A < A r < 1600 A, so a predicted continuum in the Lyman 
a forest region can then be obtained by extrapolating the fit to A r < 1216 A. 

This procedure is repeated for two sets of quasar templates. The first was generated by 
[32] using HST ultraviolet spectra of 50 low-redshift quasars, where the Lya forest continuum 
is easily estimated because of the low mean absorption. The second was generated by [33] 
from 78 SDSS DR7 ([34]) quasars with high signal-to- noise ratio, where the continua were 
fitted with a low-order spline function. The use of these different templates ensures full 
coverage of the luminosity range spanned by the BOSS quasars. The best of the two fits 
obtained from these two template sets is then selected, which we designate as Cf CA . For 
the DR9 quasars we use, the set of templates from [32] turns out to provide the best fit for 
~ 85% of the quasars. 

An additional step was applied in [31], referred to as mean flux regulation, where each 
quasar continuum was multiplied by a second order polynomial with two parameters that were 
fitted to match the mean flux evolution from [35]. This mean flux regulation substantially 
reduces the variance of the Lya absorption field, owing to the removal of spectrophotometric 
errors and of any quasar intrinsic spectral diversity that is not accounted for by the PCA 
templates. However, large-scale power in the Lya forest is also suppressed by this process 
in a way that is complex and difficult to model. For this reason, we do not apply the mean 
flux regulation procedure here. Instead, a more simple mean transmission correction (MTC) 
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is applied that is described below in section 3.3, which affects the measured cross-correlation 
with DLAs in a way that is easier to correct for. 

3.2 Mean transmitted fraction 

For the purpose of measuring the cross-correlation of the Lya forest with DLAs, it is impor- 
tant to ensure that the average value of the transmission perturbation Sp at each redshift is 
precisely zero. If this condition is not imposed, the value of the cross-correlation will differ 
from zero in the limit of large scales simply because an incorrect value of F is being used. 

For this reason, we measure the mean transmitted fraction in 50 redshift bins of Az = 
0.03 linearly spaced in redshift between z = 1.9 and z = 3.4, by computing the weighted 
average of the flux in all the pixels in each redshift bin k, centered at redshift zt, 

pPCA {z , k) = E^h/cf CA (3 2) 

where the sums are done over all the pixels i with a redshift Z{ within the redshift bin k. The 
weights are equal to the total inverse variance of fi/Cf c , 

-i 

(3.3) 

where a F = {Sp) is the intrinsic variance of the Lya forest, and F e (z) is an externally 
determined value of the mean transmitted fraction as a function of redshift, for which we use 
the result of [35]. Because these weights do not need to be obtained to very high accuracy, 
we do not iterate equations (3.2) and (3.3) to calculate the weights using the values of F PCA 
instead of F e . For the intrinsic variance as a function of redshift, we use the simple expression 

ap(z) = 0.065 [(1 + z)/3.25] 3 ' 8 , (3.4) 

based on the redshift evolution of the power spectrum measured in [36]. 

Figure 3 shows the mean transmitted fraction F PCA (red dashed line) compared to 
the result of [35] (blue dotted line). Our measurement is very close to that of [35], and is 
typically higher by only ~ 1%. In addition, there are sharp features in the inferred F{z): 
the bump at z = 2.6, and several other rapid fluctuations at lower redshift. These features 
are a systematic error that arises from the calibration of the spectra using F stars which are 
particularly important in the Balmer stellar absorption lines. Each observed plate in BOSS 
has a set of calibrations stars that are used to translate photon counts in the CCD to flux 
units [10]. The spectra of these stars are masked near the Balmer lines, and this introduces 
an artifact in the calibration vector. A similar artifact was present in earlier SDSS data 
([37]), but the effect seems to be larger in BOSS, as noted by [38]. This systematic is being 
studied and will be corrected in the future. 

The small difference between our measured F(z) and that in [35] may also be caused by 
errors in the zero-flux level computed by interpolating from neighboring fibers used for sky 
subtraction ([21]), or to other calibration systematic effects ([16]), but these errors are not a 
problem for our study of the cross-correlation as long as they are not correlated in any way 
with the presence of DLAs on nearby lines of sight. The black solid line in Figure 3 is F{z) 
in our redshift bins after the MTC is applied, which we describe next. 



F e ( Zi ) 2 ap{ Zi ) + 



(Nj) 2 
{C PCA f 
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Evolution of the measured mean flux 




Figure 3. Measured mean transmission as a function of redshift, with (solid black line) and without 
(dashed red line) the Mean Transmission Correction, compared to the measurement of [35] (blue 
dotted line). 



3.3 Mean transmission correction 

The same errors that can produce a difference in the mean transmission F PCA and the value 
F e from [35] imply the presence of stochastic errors in any single quasar spectrum that can 
systematically bias our estimated value of 8p on large scales along the line of sight. This 
may substantially increase the noise in our cross-correlation measurement. As a test of this 
effect, we measure the variance of the weighted average of 5f over the whole Lya spectrum 
of a single quasar, which we find to be on average equal to 0.03. 

The expected variance from the Lya forest alone depends on the quasar redshift, both 
because of the different comoving length of the forest and because of the evolution of the 
line of sight power with redshift. However, for z < 2.8 the expected variance is always below 
0.001. Even though metal lines and Lyman limit systems should increase this variance, 
they cannot account for a value as large as 0.03, which we infer mostly reflects the random 
variation from quasar to quasar of the systematic miscalibration. A variance of 0.03 on the 
mean of 8p over a spectrum can be induced by an rms variation of the continuum amplitude 
of ~ 17%. 

With the aim of eliminating this additional source of noise, we apply a Mean Trans- 
mission Correction (MTC) to the PCA continuum fit, which we compute for each quasar 
spectrum as 



z ieq <m [cf C A F e ( Zi )] 



where the sums are performed over all pixels i in the Lya forest region (as defined in section 
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2) of a quasar q, and the weights are computed with equation (3.3). The corrected quasar 
continuum is 



d = C 



PCA 



-4, 



9 J 



and the corrected transmission fraction at each pixel is 



F, = 



A 

Ci CF CA A n 



f, 



(3.6) 



(3.7) 



The corrected average mean transmission F(z' k ), computed as in equation (3.2) but with 
the corrected continua, is now even closer to F e (z), as shown in Figure 3. The correction 
also reduces the variance of Sp, even at the pixel level. This effect is shown in Figure 4, 
where the variance before the MTC correction is shown as the dashed red line, and after 
the MTC correction as the black solid line. When the contribution to the variance from the 
noise provided by the BOSS pipeline, after we apply a wavelength-dependent correction as 
explained in [30], is subtracted from this MTC-corrected variance, the short-dash green line 
is obtained. This line should correspond to the intrinsic variance of the Lya forest. The 
value of the variance we obtain in this way is still higher than the expected intrinsic variance 
from equation (3.4), shown as the dotted blue line, but the difference is small and is probably 
explained by residual errors in the continuum fit, metal lines and Lyman limit systems. 



A 

CM 

to 



0.1 - 



Evolution of the total pixel variance 

1 1 1 1 r 

with MTC — 
without MTC - 
measured intrinsic variance 
expected intrinsic variance 




2.2 2.4 2.6 2.8 

z 



3.2 3.4 



Figure 4. Measured variance per pixel with MTC (solid black line) and without (long-dashed red 
line). The inferred intrinsic variance using the MTC assuming the corrected noise variance provided 
by the pipeline is shown by the short-dashed green line, and the intrinsic variance assumed for the 
weights [cq. (3.4)] is the dotted blue line. 



Our estimator for the fluctuation in the transmitted flux fraction, 6f, is obtained at 
every pixel as 

Si 



CiF(zi) 



1 



(3.8) 
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where F{zi) is the mean transmission that has been obtained in the redshift bin which 
contains the pixel redshift z%. By construction, the average of 5f% at a given redshift bin is 
zero. Because F(zi) is close to F e (zi), the average of 5f% over each individual quasar spectrum 
is also nearly zero. The MTC is therefore approximately the same operation as forcing the 
weighted average of 5f to be zero on each individual quasar spectrum by subtracting its 
mean value. 

Subtracting the mean value of Sf in each line of sight removes some of the large-scale 
power of the Lya transmission field. This subtraction was also applied in [13] for measuring 
the Lya transmission autocorrelation, and an analytical expression to correct for its effect 
on the correlation function was presented in their Appendix A. We compute an equivalent 
expression of this correction to the measured cross-correlation of DLAs and the Lya forest 
in our Appendix A. 

Unless otherwise specified, our results will incorporate the MTC, and the DLA bias 
parameter will be fitted including the MTC correction of our theoretical models that is 
derived in Appendix A. However, results without including the MTC will also be presented 
and discussed. 

3.4 Estimator and covariance matrix of the cross-correlation 

The simplest method is used in this paper to obtain an estimator for the cross-correlation 
of DLAs and the Lya forest, £a'- we compute the weighted average of 5f over all the pixels 
located at a separation from a DLA that is within a bin A of r^, 



The summation symbols in this equation indicate a double sum: first over all the DLAs, and 
then over all the pixels at a separation within the A bin from the DLA. Consequently, the 
value of Wi 5 f% of one pixel appears repeated several times in the sum whenever the pixel is 
within the separation bin A from several DLAs. The DLAs are all weighted equally, and the 
Lya pixels are weighted as the inverse of their total variance, including the noise contribution 
and the intrinsic variance: 



The method is similar to the one used by [13] for measuring the Lya forest transmission 
autocorrelation. There are a number of assumptions and simplifications involved in the use 
of the simple equation (3.9) for the estimator of the cross-correlation: 

1. The correlations among the values of 5f in neighboring pixels (both along a single line 
of sight and along nearby, parallel lines of sight) can be neglected for the purpose of 
optimizing the weights. We take into account these correlations below for the purpose 
of computing the covariance matrix of the cross-correlation, but we neglect them for 
obtaining an optimal estimator by choosing the simple weights in equation (3.10). 



(3.9) 




-i -1 



(3.10) 
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2. We assume that the DLAs can be detected with a probability that is independent 
of their large-scale environment, and therefore independent of the mean Lya forest 
absorption that is superposed on their damped wings. 

Neither of these two assumptions is precisely correct. Failure of the first one makes our 
measurement suboptimal, but does not bias it in any way. The optimal manner to estimate 
the cross-correlation is discussed in Appendix B, where the assumptions involved in obtaining 
the simplified estimator in equation (3.9) are analyzed and discussed in more detail. The 
second assumption is required to avoid a systematic bias of this estimator: if DLAs are 
more likely to be detected when the forest absorption is higher, then the cross-correlation is 
artificially enhanced. We believe this second effect is small enough to be neglected, because 
tests have shown that the cross-correlation amplitude is not affected above the continuum- 
to-noise ratio that we impose on our sample (section 2.1). We nevertheless plan to make a 
correction for this effect in the future by testing how the probability of DLA detection with 
a specific algorithm is modified by the superposition of the Lya forest on the damped wings, 
using mock Lya spectra that include DLAs with a biased distribution relative to the Lya 
forest [39]. 

The measurement of the cross-correlation with equation (3.9) is done using 16 bins in 
the radial separation tt, and 8 bins in the transverse separation a. To test for the symmetry of 
the cross-correlation under a sign change of tt, we use different bins for negative and positive 
7r, where tt is positive when the pixel of the Lya forest transmission perturbation is at higher 
redshift than the DLA. The bins in tt are, in units of comoving /i" 1 Mpc, bounded by the 
values (—60, —40, —30, —20, —15, —10, —6, —3, 0) and the same positive values, while 
the bins in a are bounded, in the same units, by (1, 4, 7, 10, 15, 20, 30, 40, 60). Pairs at 
transverse separation a < 1 h~ l Mpc are not used (the number of these pairs is negligibly 
small in our sample). This results in a total of 128 separation bins in a and tt. The weighted 
average values of (tt, a) of all the contributing pixel-DLA pairs to every bin are stored, with 
the same weights as in equation (3.10 (these are generally close but not exactly equal to the 
central values of each bin). The measurement is generally done using a single bin of the 
mean redshift z of the Lya forest pixel and the DLA, which is required to be in the range 
2.0 < z < 3.5, although we also present some results in the next section using three redshift 
bins. 

The covariance of the cross-correlation measured with equation (3.9) in two bins A and 
B is equal to 



where Cy = (^Fi^Fj) is the correlation of the measured values of 5f in pixels i and j, 
separated in redshift space by rjj. Note that the correlation Cij, with two subindexes for 
the two correlated pixels, should not be confused with the quasar continuum which always 
has one subindex, and the covariance matrix of the cross-correlation, C, is a different matrix 
with indexes referring to bins in (ir, a). There are three main contributions to the correlation 
Cij: the intrinsic autocorrelation of the Lya forest at a given separation ^(r^), the noise 
term (-/V?) / (CfF?) that we assume to be uncorrelated among different pixels, and continuum 




(3.11) 
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fitting errors with a correlation £ c ( r ij) that affect only pairs of pixels in the same spectrum, 



Here, is the perpendicular component of r^, and the Dirac delta function 8 indicates 
that the last term is non-zero for pixels on the same quasar spectrum only, as opposed to the 
Kronecker function bf- in the second term that indicates that the noise term is only present 
when both pixels are the same. 

In general, a joint analysis of the Lya auto-correlation and a cross-correlation may be 
done, and the measured correlation £p can then be used to compute the cross-correlation 
covariance. Here, we use instead a theoretical correlation function computed as the Fourier 
transform of the Lya power spectrum measured in [40] from numerical simulations, with 
modified bias parameters in agreement with the measurement of the correlation function 
measured with the first year of BOSS data in [13]. 

We use a Lya density bias parameter 6p = —0.168 and a redshift distortion parameter 
/3f = 1 at z = 2.25, and we include the non-linear term D(k,Hk) with the parameter values 
of the fiducial model in the first row of Table 2 of [40]. These values are consistent with 
the results of [13] . The amplitude of the correlation function is assumed to be proportional 
to (1 + z) a , with a = 3.8, as found from the ID power spectrum measurement in [36]. We 
expect most of the correlated errors in the continuum fitting to be removed by the MTC 
procedure, so we generally do not include the term £ c in equation (3.12) except when the 
MTC is omitted. In the latter case, a constant contribution £ c = 0.03, equal to the measured 
variance of the value of 8p averaged over one spectrum, is added to C%j for all pixel pairs on 
a common line of sight. 

The calculation of the cross-correlation covariance using equation (3.11) can be time 
consuming because of the need to evaluate £p for all pixel pairs with separations within the 
A and B bins from every pair of DLAs. Each sum in the numerator of equation (3.11) is 
over every DLA, and over every pixel within the separation bin. To increase the speed of the 
calculation, we include only pixel pairs with a perpendicular separation smaller than oij = 
20 h^ 1 Mpc, a similar approximation to the one used in [13] for computing the covariance 
matrix. This approach ensures the inclusion of the most important correlations of our Lya 
spectra in the covariance matrix. 

3.5 Sub-samples and bootstrap errors 

The errors obtained on the measured cross-correlation using the covariance matrix computed 
with equations (3.11) and (3.12) rely on the accuracy of this calculation, which may be 
affected by the assumed Lya autocorrelation and the neglect of continuum fitting errors 
and other possible systematics (such as correlated sky subtraction errors). It is therefore 
useful to compute an alternative set of bootstrap errors by dividing our quasar sample into 
a number of sub-samples and testing the variation of the obtained cross-correlations among 
the sub-samples. 

We divide the area covered by DR9 into 12 sub-samples as shown in Figure 5. Table 1 
provides the number of quasars and the number of DLAs present in each sub-sample. The 




(3.12) 
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DLA-Lya cross-correlation and its covariance matrix are separately computed in each sub- 
sample. The typical sub-sample size is ~ 10 deg, or ~ 700 h~ l Mpc , which is much larger than 
the maximum scale of our analysis, a < 60/i _1 Mpc. We are therefore not concerned about 
the loss of a small fraction of DLA-Lya pairs that are near the border of two neighboring sub- 
samples. Because of the geometry of the observed sky area in DR9, most of the sub-samples 
have short borders with their neighbors. 



Definition of the 9 northern sub-samples Definition of the 3 southern sub-sampies 




100 150 200 250 300 300 320 340 360 380 400 420 

J2000 RA (deg) J2000 RA (deg) 



Figure 5. The 12 sub-samples in which we split the DR9 quasar sample, 9 in the northern Galactic 
hemisphere (left) and 3 in the southern one (right). 



Sub-sample 


# quasars 


# DLA 


N-0 


2290 


303 


N-l 


5176 


647 


N-2 


5264 


670 


N-3 


5246 


758 


N-4 


4010 


537 


N-5 


4056 


545 


N-6 


4823 


658 


N-7 


3255 


406 


N-8 


6322 


819 


S-9 


4339 


543 


S-10 


4308 


503 


S-ll 


3360 


391 



Table 1. Number of quasars and DLAs included in our analysis in each of the 12 sub-samples. 

The cross-correlations and their covariance matrices are obtained in each sub-sample, 
£ Q and C a , and are combined using the equations 

C^E^ 1 ' ( 3 - 13 ) 

a 

i=CY / C~ 1 L- (3.14) 

a 

Bootstrap errors are calculated from these 12 sub-samples by generating N = 100 
random combinations of the 12 sub-samples with repetitions, and recomputing C and £ for 
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each of these 100 combinations [41]. The dispersion among these combinations of the values 
of any model parameter that is fitted to the result yields the bootstrap error. 



3.6 Fitting the DLA bias 

In the limit of large-scales, linear theory predicts the precise form of all the correlations 
among any pair of tracers of the large-scale structure in redshift space. The Lya forest probes 
material at low density and high redshift, which is less affected by non- linear evolution and 
random peculiar velocities than a set of galaxies moving on orbits within gravitationally 
collapsed halos. Therefore, linear theory remains a good approximation down to smaller 
scales for the Lya forest than for other tracers. In linear theory, redshift distortions cause 
the amplitude of each Fourier mode in a biased tracer field to be enhanced by the factor 
6(1 + /3 nl), where b is the density bias factor, (3 is the redshift distortion parameter, and 
is the cosine of the angle between the Fourier mode and the line of sight [42]. The linear 
cross-power spectrum of the DLAs and the Lya forest is therefore equal to 

P DF (k, z) = b D (z) [1 + p D {z)iJ%\ b F {z) [1 + Pf(z)hI] P L {k, z) . (3.15) 

where bo and b F are the DLA and Lya density bias factors, f3o and f3 F are their redshift 
space distortion parameters, and P F {k,z) is the linear matter power spectrum. The cross- 
correlation function is the Fourier transform of Pdf, an d can be computed using the following 
equations (derived as in [43]) 

£(r) = £o(r)P (/i) + 6(r)P 2 (/i) + U{r)P^) , (3.16) 

where fi = ir/r is the angle cosine in redshift space, Po, P2 and P4 are the Legendre polyno- 
mials, and the functions £oj £2 and £4 are 

&(r) = b D b F [1 + ({3d + 0*0/3 + Pd{3 F /5] C(r) , (3.17) 

6(r) = b D b F [2/3((3 D + (3 F ) + A/7(3 D [3 F ] [C(r) - C(r)] , (3.18) 

U(r) = 8/35 b D b F p D p F [C(r) + 5/2C(r) - 7/2<f(r)] . (3.19) 

The function £(r) is the linear correlation function in real space, and and C,(r) are 
averages of £(r) within r given in equation (9) of [43]. We compute £(r) using the ACDM 
model with the parameters given at the end of the Introduction. 

We assume the DLA bias 6d to be constant with redshift, while for the Lya forest, we 
assume the bias evolves as b F (z)/b F (z = 2.25) = [(l + z)/3.25] a//2+1 , where a = 3.8 describes 
the evolution of the amplitude of the line-of-sight power spectrum measured in [36]. The 
rapid evolution of the amplitude of 5 F in the Lya forest is a consequence of the rapid change 
in the mean transmitted fraction with redshift. 

For the redshift distortion of the DLAs, we assume (3d = f(p,)/bo ~ b^ 1 , where /(O) is 
the logarithmic derivative of the linear growth factor [42], which is close to unity at z > 2 in 
the ACDM model. This assumption may not be exactly correct for DLAs if the probability 
of detection of a DLA depends on the large-scale peculiar velocity gradient because of the 
Lya forest absorption that is superposed on the damped wings. This effect would generate 
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a peculiar velocity bias analogously to the case of the Lya forest, but we expect this effect 
to be small. 

For the Lya forest, bp and f3p are independent because the 5f tracer is subject to a 
strong bias by the large-scale peculiar velocity gradient [40]. We use the results for the Lya 
autocorrelation of [13] to fix 6f(1 + /3p) = —0.336. Most of our results will be presented 
assuming (3p = 1-0, in which case the only free parameter is bo- Some results will also be 
obtained when /3p, which was poorly determined by [13], is also treated as a free parameter. 
Note that the Lya forest has a negative bias, since denser regions of the universe have a lower 
transmitted flux. This implies that the cross-correlation is also negative. 

The model of the linear theory in equation (3.16) is generally modified to account for the 
effect of the MTC, which removes the mean value of Sf in each quasar spectrum, therefore 
distorting the DLA-Lya cross-correlation. The way this correction is modeled and calculated 
is explained in Appendix A. Instead of evaluating the theoretical model at the center of 
every bin in the a, tt components of the DLA-Lya separation, we evaluate it at the weighted 
average values of a, tt for all the DLA-Lya pairs that contribute to the bin. To minimize 
possible non-linear effects on the correlation function, all the fits are done using only bins 
with r > 5 h^ 1 Mpc; this eliminates the two bins with smallest a and tt, reducing the number 
of bins to 126. 

For each theoretical model that predicts the cross-correlation at each bin A of (vr,cr), 
£,t,A, we compute its \ 2 as 

x 2 = Ha - itA^Aliis - e t ,B) , (3.20) 

where the indexes A and B are summed over all bins. The likelihood function is 

exp[-0.5v 2 ] 

L= \ T/ - J , 3.21 
{2tt) n / 2 \C\ K ' 

where |C| is the determinant of the covariance matrix of the cross-correlation (computed 
using equations 3.11 and 3.13), N = 126 is the number of bins, and we have assumed that 
the distribution of errors in the cross-correlation is Gaussian. 

Our best estimate of the DLA bias is the one that maximizes the likelihood, and its 
uncertainty is evaluated using two different approaches. First we estimate the errors on the 
fitted parameters with the Monte Carlo Markov Chain technique (MCMC errors), where 
we sample the parameter space using the likelihood function described above and compute 
the dispersion around the best fit value. Secondly, the fits are repeated with the bootstrap 
realizations of the 12 sub-samples randomly selected with repetitions as described previously, 
obtaining the bootstrap error (BS error) from the dispersion in the fitted parameters among 
the bootstrap realizations. Both the MCMC and BS errors will be presented. 

4 Results 

In this section the results of the cross-correlation of DLAs and the Lya absorption in redshift 
space are presented. After obtaining a fit to our fiducial model with only one free parameter, 
the bias factor of DLAs, we analyze the consistency of the covariance matrix and the bootstrap 
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errors, the evidence for the presence of redshift distortions, and the dependence of our best 
fit model on the range of separations that are used. We then investigate a variety of possible 
dependences of the DLA bias factor on several variables: redshift, column density, and the 
Lya forest redshift distortion parameter. In general the results presented here have been 
corrected with the MTC, and the fitted models are also corrected as described in Appendix 
A, except when we test for the effect of the MTC. 

4.1 Measured cross-correlation and bias parameter 

The results of the cross-correlation as a function of tt for our fiducial model are shown as 
red points with errorbars in figures 6 and 7 in eight panels, corresponding to our eight bins 
in a. The errorbars are the square root of the diagonal elements of the covariance matrix 
C computed as described in section 3.4 (equation 3.13). The blue solid curves are our best 
fit fiducial model (all of our fits exclude bins at r = (a 2 + 7T 2 ) 1 / 2 < 5/i _1 Mpc, which are 
only the two bins at the lowest a and tt, even though they are shown in the figures), which 
fixes (3f = 1 and fits only one free parameter, the DLA bias. The result obtained from our 
covariance matrix is = 2.17 ± 0.20, with a x 2 of 106 for 125 degrees of freedom, and the 
bootstrap analysis gives the same errorbar. Note that we are measuring only the amplitude 
of the cross-correlation, which is proportional to &d&f<t| and depends on /3p, and the error on 
br> here reflects only that in our measurement of the amplitude. The results are also shown 
in a contour plot in the left panel of figure 8, and the best fit fiducial model is shown in the 
right panel with the same contour levels. 

As seen in the figures and indicated by the x 2 value, the fit of the standard ACDM model 
to the DLA-Lya forest cross-correlation is excellent. The measurements are fully consistent 
with the radial dependence and anisotropy that is expected in linear theory for the model 
that agrees with all other cosmological determinations of the large-scale power spectrum. 

4.2 Effect of the MTC 

All our results are generally obtained including the MTC, and the fitted models include the 
correction described in Appendix A. We now check the effect of the MTC by recomputing 
the cross-correlation without including the correction, i.e., by using directly the PCA-only 
continuum. In this case, we add to our covariance matrix the additional term £ c = 0.03 in 
equation (3.12) due to the error of the continuum fitting for pixels on a common line of sight. 

The measurements with and without the MTC are compared in figure 9 (red and black 
errorbars), for four selected a bins. The best fit model is also shown with and without 
including the correction of appendix A (blue solid and dotted green line, respectively). Table 
2 shows the result for the bias for the NOMTC case (uncorrected data and uncorrected fit 
model), and for the case where we apply the MTC to the data but we leave the fitted model 
uncorrected (NOCOR case). 

The effect of the MTC is clearly a minor one; it has little effect on the result on the 
DLA bias, and the errorbars are only increased by ~ 10% when the continuum fitting errors 
are included. The effect of the correction of appendix A is small compared to the present 
errorbars of the measurements, but our model characterization of the effect of the MTC is 
supported by the larger value of % 2 in the NOCOR case. We keep the correction despite 
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Figure 6. Measured cross-correlation in different bins of a: from top-left to bottom-right, 1 < a < 
4 h^ 1 Mpc. 4 < a < 7 h^ 1 Mpc, 7 < a < 10 ft -1 Mpc, and 10 < a < 15 h^ 1 Mpc. Solid lines show the 
best fit model for the fiducial analysis, including the correction derived in appendix A. 





b D 


BS errors 


MCMC errors 


X 2 (d.o.f) 


FIDUCIAL 


2.17 


0.20 


0.20 


106 (125) 


NOMTC 


2.11 


0.21 


0.22 


109 (125) 


NOCOR 


2.00 


0.19 


0.20 


111 (125) 


NODLA 


2.25 


0.22 


0.21 


109 (125) 


LOWZ (2 < z < 2.25) 


2.18 


0.41 


0.33 


116 (125) 


MIDZ (2.25 < z < 2.5) 


2.16 


0.32 


0.34 


109 (125) 


HIGHZ (2.5 < z < 3.5) 


1.88 


0.57 


0.37 


92 (125) 


LOWNHI (log(jVtfj) < 20.4) 


2.27 


0.30 


0.29 


133 (125) 


HIGHNHI (log(N HI ) > 20.4) 


1.89 


0.26 


0.30 


110 (125) 



Table 2. Best fit value of the DLA bias with its bootstrap errors, MCMC errors and x 2 value of the 
fit, for various analyses: FIDUCIAL (with the MTC and the theory corrected with the expression 
derived in appendix A), NOMTC (PCA-only continuum fitting, uncorrected theory), NOCOR (MTC, 
uncorrected theory), NODLA (spectra containing DLAs are rejected), data split in redshift bins 
(LOWZ, MIDZ, HIGHZ) and finally the DLA sample split in two bins of column density (LOWNHI, 
HIGHNHI). 
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Figure 7. Same as figure 6, for the a bins, from top-left to bottom-right, 15 < a < 20 h 1 Mpc, 
20 < a < 30 hr 1 Mpc, 30 < a < 40 hr 1 Mpc, and 40 < a < 60 hr 1 Mpc. 
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Figure 8. Contour plot of the measured redshift space cross-correlation for our fiducial analysis (left) 
and our best fit theoretical model that includes the MTC correction (right). 



the fact that with the present level of noise its effect is not significant, because it will be 
important for future studies with smaller errorbars and larger separations. 
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Figure 9. Effect of the MTC on the measured cross-correlation in the four selected bins of a that are 
indicated. The fiducial measurement is plotted with red errorbars, while the measurement without 
the MTC is shown in black (with a small horizontal shift for better visualization). The blue solid 
line is the fiducial theory including the MTC from appendix A, and the dotted green line without the 
MTC. 



4.3 Tests of the covariance matrix 

We have seen that the bootstrap error on the fitted bias parameter of the DLAs is generally 
in good agreement with the one computed from our covariance matrix using the MCMC 
technique (see table 2). As a further test of the accuracy of our errors, we now compare the 
scatter between the values of the measured cross-correlation in bins, using the 12 sub-samples 
defined in section 3.5, with the predicted errors from the diagonal elements of the covariance 
matrix. 

For the purpose of this comparison, instead of averaging the cross-correlation values 
from the sub-samples weighted with their inverse covariance matrices, as in equation 3.14, 
we average by weighting each sub-sample with the sum of the weights of all the pixels that 
contributed to each bin A, w^, where the index a labels each sub-sample. This estimate of 
the cross-correlation in bin A is 

U = "M" ' (41) 
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and the uncertainty a a of this estimate is 
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Figure 10. Errors of the cross-correlation in four selected bins of cr, as indicated in each panel. The 
solid lines show the estimate from the scatter between sub-samples (equation 4.2), while the dashed 
lines show the estimate from the diagonal elements of the covariance matrix. 

In figure 10 the errors computed from the scattering among the sub-samples and from 
the diagonal elements of the covariance matrix are compared. The errors from the scatter 
among sub-samples (solid lines) are noisy because of the small number of sub-samples, but 
in general they agree remarkably well with the ones obtained from the covariance matrix 
(dashed lines). The thick lines are for the fiducial case, and the thin lines for the NOMTC 
case. As before, the NOMTC case incorporates continuum fitting errors that account for the 
larger covariance errors by ~ 10%. Generally the NOMTC case has larger scatter than the 
MTC case, as we would expect if the MTC actually helps reducing the uncertainty in the 
cross-correlation measurement. 



4.4 Evidence for the redshift distortions 



The gravitational evolution origin of the DLA-Lya absorption cross-correlation we detect is 
tested by the presence of the predicted redshift distortions in linear theory. The expected 



elongation of contours perpendicular to the line of sight is apparent in figure 8. By fitting only 
the DLA bias factor and making the reasonable choice @f = 1 based on previous observations 
of the Lya autocorrelation in [13], the observed anisotropy is consistently matched by our 
linear model. The detection of redshift distortions can be quantified by comparing with a fit 
that forces isotropy requiring (3p = /3d = 0. This isotropic fit yields x 2 = 150, implying a 
detection of more than 6a when compared to \ 2 = 106 (125 d.o.f.) for our best fit model with 
/3f = 1 and (3d = /(^)/&D (points at r < 5/j" 1 Mpc are not used in the fits in both cases). 
This result provides strong support for the interpretation that the measured cross-correlation 
is induced by the gravitational growth of large-scale structure that is traced by the DLAs 
and the Lya forest. 

4.5 Scale dependence of the DLA bias 

In addition to the redshift distortions, the predicted dependence of the cross-correlation on 
the separation r in the ACDM model is another important test that our measurement agrees 
with the theoretical expectation. This can be rephrased as testing that the inferred DLA 
bias factor is constant in different intervals of r. The results of this exercise are shown in 
figure 11. 
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Figure 11. Fitted DLA bias in several bins of the separation r for our fiducial analysis. The green 
lines show the values obtained when combining all bins above r > 5 h^ 1 Mpc. 

The value of the bias is consistent with being constant for r > 5 h^ 1 Mpc. At smaller 
scales we expect non-linear effects to be important, so it is not surprising that the first point 
has a discrepant value of the bias. The first bin that we use in our fits, 5/i _1 Mpc < r < 
7/i _1 Mpc, also has a lower bias than the ones at larger radii by 1.3er, and it may already be 
affected by non-linear effects to some degree. In addition to the physical non-linearities in 
the distribution of gas and DLAs in redshift space, the transformation from optical depth to 
the transmitted fraction introduces an additional non-linearity that in this case may be the 
dominant effect. 
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The presence of non-linearities may be accounted for in the future by using predic- 
tions from hydrodynamic numerical simulations of structure formation. When the full BOSS 
dataset is available, the reduced error bars will allow for a better measurement of these 
non-linear effects. We have also checked the symmetry of the cross-correlation under a sign 
change of the tt coordinate. There is no statistically significant difference between the fitted 
fiducial model when only the bins with tt > or those with tt < are used. 

4.6 Effect of the DLAs on the Lya forest 

In our fiducial analysis we include in the Lya sample 3047 lines of sight that contain at least 
one DLA. As explained in section 2, we apply a mask on the central region of the absorption 
profiles and correct for the damped wings by multiplying the continuum estimate with a 
Voigt profile. To test the effect of this correction we measure the cross-correlation using 
only the remaining 49402 lines of sight. The results can be seen in table 2 under the name 
NODLA, and are clearly consistent with our fiducial analysis. 

4.7 Dependence on the Lya bias 

Our fits for the DLA bias are based on assuming fixed values for the Lya forest bias param- 
eters: the well constrained combination 6f(1 + (3f) = —0.336 from [13], and (3p = L The 
value of Pf is still highly uncertain, but improved measurements are expected in the near 
future from BOSS. We therefore examine the dependence of our result for bn on the assumed 
value of /3f- 
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Figure 12. Best fit value of the DLA bias for different values of f3p, for our fiducial analysis. The 
dependence is well described with the fitting formula bo = (2.17 ± 0.2)(3p 22 (blue line). 

In figure 12 we plot the results assuming different values of j3p. We show that for the 
most plausible range 0.4 < (3p < 2.0, the dependence is well fitted by a simple function 

b D = (2.17 ±0.20) • ( 4 -3) 
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Figure 13 shows the likelihood contours for the two-parameter fit when 6p is allowed 
to vary, keeping fixed the combination 6f(1 + /3p) = —0.336, for our fiducial analysis. The 
data clearly prefer larger values of (5f (smaller values of bp) than the one we have assumed 
so far (/3f = 1, bp = —0.168). However, the best fit value of @f is quite sensitive to the 
applied correction, and so we prefer to be conservative by assuming a fixed, low value of 0f, 
minimizing the required value of the DLA bias. 
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Figure 13. Degeneracy between the value of bo and 6p when fitting both parameters at the same 
time in our fiducial analysis. Since we keep fixed the well constrained quantity 6^(1 + ftp) = —0.336, 
each value of bp is related to a value of ftp. The horizontal green line shows the best fit value when 
fixing j3p = 1 (frp = —0.168), bo = 2.17. The blue line shows the fitted function in equation 4.3. 

4.8 Redshift evolution 

Until now we have assumed that the value of the DLA bias is constant with redshift. We 
now evaluate the cross-correlation in three different redshift bins, 2.0 < z < 2.25 (LOWZ), 
2.25 < z < 2.5 (MIDZ) and 2.5 < z < 3.5 (HIGHZ). The third bin is wider because the Lya 
forest data available for cross-correlation is much more sparse at high redshift in BOSS. We 
assume the redshift evolution of the Lya forest bias factor to be ftp oc (1 + z) 2 ' 9 , based on 
the evolution of the line of sight power spectrum measured in [36]. The best fit values for the 
three redshift bins, shown in table 2, provide some evidence that the DLA bias decreases with 
redshift. This evolution is opposite to that expected if the DLA host halos are approximately 
at constant velocity dispersion: in this case the DLA bias should increase with redshift as 
the host halos become increasingly rare fluctuations. The observed trend is, however, of low 
statistical significance. 
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4.9 Dependence on column density 

Finally, we split our sample of DLAs into two data sets as a function of column density, 
separated at Njji = 10 20 ' 4 cm -2 . If systems with higher column density resided in halos of 
higher mass, one would expect the bias of the high column density sub-sample (HIGHNHI) 
to be larger than that of the low column density (LOWNHI). The results shown in table 2 
are consistent with a bias that is independent of column density, showing, if anything, an 
opposite trend that is not statistically significant. 

A possible systematic error that may affect our DLA bias is the inclusion in our DLA 
catalogue of systems that are not actual DLAs, but regions of absorption with lower column 
density than our threshold which pass our column density cut because of a combination of 
spectral noise and the Lya forest superposed on the damped wings. At the same time, DLAs 
that are just above our column density cut may more easily be attributed a column density 
below the threshold (and be therefore omitted from the catalogue) when the superposed 
Lya forest is particularly weak. If this error is important, the Lya forest superposed on the 
wings of DLAs may introduce a correlation with the Lya forest in nearby lines of sight that 
would systematically enhance our measured bo- We do not believe this effect is important 
because preliminary tests of the DLA selection with mocks show that there is a low rate 
of miss-identifications ([26]), and because of the consistency of the radial dependence of 
the cross-correlation and the anisotropy with theoretical expectations found in sections 4.4 
and 4.5, but we plan to examine this selection effect more carefully in the future using the 
techniques described in [39]. 

4.10 Result for the sample affected by BAL Contamination 

As seen in figure 1, there is a significant excess of systems detected in the window 1005 A < 
A r < 1037 A, probably due to contamination by BALs with small Balnicity index BI. The 
systems detected in this window were rejected in the fiducial analysis. We measured the 
cross-correlation using only the 989 systems that fall into this range of restframe wavelength, 
including those detected in lines of sight with a Balnicity index in the range < BI < 
1000 kms" 1 . The measured bias is bo = 0.68±0.50 (BS errors), showing that at least half of 
these systems are indeed contaminants, and confirming that the measured cross-correlation 
is present only for physically real DLAs. 

5 Discussion 

The measurement of the cross-correlation of DLAs and the Lya forest absorption, presented 
for the first time in this paper, provides a new observational constraint on the nature of 
DLAs: the bias factor of their host halos. We have measured this bias factor to be bo = 
(2.17 ± 0.20) (3p 22 , at a mean redshift of z = 2.3. The analysis of numerical simulations of 
the Lya forest predict a value of /3p ~ 1.5 ([40]), and even though the first measurement of 
/3p from the Lya forest autocorrelation [13] allowed for a wide range of /3p around unity, a 
more recent analysis also suggests a value larger than one (Slosar, private communication). 

As discussed in the introduction, the bias factor constrains the mass distribution of the 
host halos of DLAs. A small bias factor implies low halo masses and an association of most 



- 25 - 



DLAs with dwarf galaxies, while a large bias factor means that most DLAs are extended 
clouds in massive halos. The observational determination of the mean bias factor is therefore 
a new constraint that models of DLAs must satisfy. 

Theoretical models of DLAs need to provide the average cross section of a halo of 
mass M for producing an absorption system with Nhi > 10 20 ' 3 cm -2 . The predicted rate 
of incidence is given by equation (1.1), and the new observable of the mean bias factor 
to be matched by models is obtained with equation (1.2). Numerical simulations of galaxy 
formation following the collapse and radiative cooling of gas in halos have found that the rate 
of incidence predicted for DLAs is affected by two uncertain factors: the limited resolution 
of the simulations and the impact of galactic winds ([44], [45], [46]). As the resolution is 
improved and halos are well resolved down to the lowest masses for which photoionized gas 
can collapse, the rate of incidence is increased; at the same time, when a prescription for 
supernova-driven galactic winds is included leading to the ejection of gas from low-mass 
halos when star formation occurs, the gas mass and the DLA cross sections in low-mass 
halos decrease as the energy that is deposited in winds is increased ([47], [48], [49], [50]). 

These models based on cosmological simulations do not at present provide fundamental 
physical predictions for the properties of DLAs, because the way that the gas distribution 
in halos is affected by cooling and fragmentation, radiative transfer, star formation and 
galactic winds depends on extremely complex physical processes that can only be roughly 
approximated through simple parameterized recipes that are incorporated into the equations 
being solved in the simulations. Nevertheless, they can provide a basic guide as to the 
general characteristics of star formation rates and wind strengths that are needed to satisfy 
the observational constraints. 

5.1 Constraints on the halo mass distribution from the bias factor 

We now examine the expected bias factors for simple relations of the DLA cross section and 
the halo mass that are useful fits to the results found in the numerical simulations of DLAs. 
The DLA systems are likely distributed over a wide range of halo mass, but they are not 
expected in halos that are too shallow to hold photoionized gas and allow it to radiatively 
cool. This condition corresponds to a halo circular velocity of v c ~ 20kms _1 , or a halo mass 
at z = 2.3 M h ~ v^/(3ttHG) ~ 10 9 M . Many numerical models of DLAs, starting with [45], 
have obtained fits to simulation results with a power-law relation for the cross section, 

S(M) = S (M/M min ) Q (M > M min ) . (5.1) 

The distribution of halo masses and the bias factor of halos has been thoroughly examined in 
analytic models and numerical simulations and are robustly predicted in the ACDM model 
(e.g., [6]; [7]). Here we use the mass distribution and halo-bias relation given by [7]. 

In the left panel of figure 14, the halo mass-bias relation of [7] is shown as the solid 
purple line. The other lines show the average value of the bias weighted according to the 
power-law cross section of equation (5.1), for several values of a, as a function of M m i n . All 
the relations are computed at z = 2.3, the mean redshift of our DLA bias measurement. Halos 
having the value of the mean bias we have measured for DLAs have masses of ~ 6 x 10 11 M & 
(we use h = 0.71). Realistically, the host halo masses should vary over a broad range, but 
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Figure 14. Average bias as a function of M m i n when the cross section depends on halo mass as M£ 
with a sharp lower limit at M m i n (left), and for the model of equation (5.2) with a smooth suppression 
at masses below Mq, for several values of a, computed at z = 2.3 using [7]. The horizontal dotted 
lines show the result of this paper for the mean DLA bias with la errorbars. 

this mass should be a typical one. Our measurement of the bias factor clearly requires DLAs 
to be mostly associated with host halos of relatively massive galaxies. 

Fixing the lower limit in equation (5.1) at M m j n ~ 10 9 Mq, the best value of a that 
fits our measured bias for / 8jr = lisa = l.l±0.1. To better interpret the implication of 
this result, consider a simple model where the DLAs are clouds found within a fixed fraction 
f a of the virial radius of the halo, r„j r oc M 1 / 3 , with a covering factor f c . Then the cross 
section of a halo is E(M) oc /^/ C M 2 / 3 . In other words, a = 2/3 if DLAs are found within 
a fixed fraction of the virial radius and have a constant covering factor. This value of a is 
ruled out by our result at the 4<r level for M m j n = 10 9 Mq. Therefore, DLAs must cover an 
increasing fraction of the projected area of halos with increasing mass, up to masses larger 
than ~ 10 12 Mq. For f3p = 0.6, a = 2/3 is still ruled out at the 3a level. 

Most of the models from cosmological simulations predict shallower slopes than are 
implied by our result of the DLA bias. Models without strong effects from galactic winds 
predict a slope close to a = 2/3 (see, e.g., [48], [50], [51], [52], [53] for recent work on 
the effects of galactic winds and differences among hydrodynamic codes). Models with the 
strongest winds are barely consistent with our measured bias factor: model Q5 of [47], [48], 
with a = 1 and M m j n ~ 10 9 Mq, yields a bias about 1.5cr below our measurement. The model 
of momentum-driven winds of [51], with a = 0.92 at z = 2.25, is also barely consistent. 

The right panel in figure 14 shows the same relation of the mean bias for a model where 
the sharp cutoff of the cross section at M m j n is replaced by a smooth function: 

We find that the numerical results of [50] (see their figure 4) are well matched by this formula 
with a = 0.5, M = 10 9 - 5 M Q , and E = 40kpc 2 . Their results are at z = 3; to shift them 
to z = 2.3 we assume that the properties of DLAs stay constant in halos of fixed circular 
velocity, which shifts their characteristic halo mass to M = 10 9 - 7 M . Their model predicts 
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a mean bias factor bo ~ 1.35, which is ruled out by our measurement at more than the 3a 
level if [ip > 1- Even if (3p = 0.6 their model is ruled out at a confidence level of 2.5a. 

Note that the models need to predict only the average DLA cross section as a function 
of the halo mass, because the mean DLA bias does not depend on the scatter of the cross 
section for a fixed halo mass. 

5.2 Constraints on the cross sections from the rate of incidence 

The most recent determination of the rate of DLA incidence using the BOSS 9th Data Release 
shows that a fraction 0.226 ± 0.005 of the sky is covered by a DLA with Nhi > 10 20 ' 3 cm -2 
in the redshift range 2 < z < 3, as derived from the column density distribution in table 1 of 
[26] (the error quoted here is only the statistical one obtained from the same table, and docs 
not include the systematic error due to corrections for incompleteness and impurity). This 
corresponds to a rate of incidence per unit of absorption pathlength, dx = (1+ z) 2 Hq / H (z) dz, 
equal to 0.066 ± 0.002, which we assume to be constant with redshift This observed rate of 
incidence fixes a required normalization of the relation of the DLA cross section and halo 
mass for any model. 

We consider here a few examples and examine the value of the physical cross section 
£ at the characteristic halo mass = 10 12 Mq to illustrate the typical implications of the 
observed rate of incidence and bias factor of DLAs: 

1. A model with a = 1 and M m j n = 1O 9 M0, which produces a bias bp = 1.8, requires 
£(M = 10 12 Mq) = 1400 kpc 2 to match the observed incidence rate. This is a factor of 
two larger than the cross section in the Q5 model of [47], and a factor 1.5 larger than 
the prediction of the momentum-driven wind model of [51]. 

2. A model with equation (5.2), with a = 1 and Mq = 10 10 Mq, matches our observed 
bias for S(M = 10 12 M Q ) = 2400 kpc 2 . 

3. A model with equation (5.2), with a = 0.5, needs an increased Mq = 3 x 10 11 Mq to 
match the observed bias, and then requires £(M = 10 12 M Q ) = 4400 kpc 2 . 

This illustrative example shows that the measured DLA bias factor can be matched either by 
an extended DLA host halo distribution with a > 1 and a low mass cutoff, or with a narrower 
range of halos around the characteristic mass of 10 12 Mq. A more extended distribution of 
halo mass reduces the required cross section because the rate of incidence is accounted for 
by a wider range of halo masses. 

5.3 Consequences for the nature of DLAs 

Generally, the models of DLAs discussed in the literature based on hydrodynamic simulations 
of structure formation fail to correctly predict the properties of DLAs for a common reason: 
the cross sections they predict for massive halos (M ~ 10 12 M Q ) are too small. The models 
that come closest to matching the observations are the Q5 model of strong winds in [47], 
which still predicts a rate of incidence and a baryonic content of DLAs that is a factor of 2 
below the observed ones, and the momentum-driven wind model of [51], which with a = 0.92 
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would predict a DLA bias that is still 2a below our best value for ftp = 1, and cross sections 
that are too small by a factor ~ 1.5. 

Our observational result therefore implies that some neglected physical mechanism keeps 
the atomic gas spread out over an area larger than predicted by present models in host halos 
of massive galaxies. The natural expectation is that galactic winds that are even stronger 
than in the Q5 model of [47] and momentum-driven wind model of [51], and that have an 
increased impact on the halo gas distribution, are required. We emphasize, however, that 
DLAs in low-mass halos do not need to be further suppressed: in fact, these objects also 
help account for the total rate of incidence. The requirement for stronger winds is to allow 
atomic gas to stay at large radius over several orbits in massive halos. 

To better understand the consequences of our result, we consider the Q5 model of [47] 
with all the cross sections increased by a factor two, to match both the bias factor and rate 
of incidence of DLAs. Then, the mean cross section in a 10 12 M & halo with circular velocity 
v c ~ 200 kins" 1 at z = 2.3 is E = 1400 kpc 2 . The corresponding proper radius containing 

— 1/2 

the DLAs for a circular cross section is ~ 20/ c kpc, or about 20% of the virial radius of 
the halo if the covering factor f c is near unity, and the orbital period is ~ 5 x 10 8 / c years. 
These clouds can therefore complete no more than a few orbits at this radius over the age of 
the universe at z = 2.3, and would be forming a halo system that is dynamically supported 
by random motions, because the large radius implies that rotation may account for only a 
small fraction of the velocity dispersion. 

The massive nature of most DLA host halos helps account for the large velocity widths 
observed for the metal lines associated with DLAs ([54], [55]), as noted also by [56]. At the 
same time, the problem of the rate of energy dissipation in DLAs ([57], [49]) is solved, because 
the long orbital periods imply long collision times among the DLA clouds, and the available 
gravitational energy that can be dissipated in cloud collisions is increased in massive halos. 

Once the cross sections are known, we can infer the fraction of baryons present in 
DLA clouds if we make the simple assumption that the distribution of column densities is 
independent of halo mass, which is supported by the lack of dependence of our measured bias 
factor on column density. From table 1 in [26], we infer that the mean column density of all 
DLAs with N m > 2 x 10 20 cm" 2 is N HI = 7.8 x 10 20 cm~ 2 . In a halo of 10 12 M of total 
mass and cross section E = 1400 kpc 2 (and using a hydrogen abundance by mass X = 0.76), 
this implies that a fraction of 8% of all the baryons in the halo are in the form of atomic 
gas in the DLA clouds. Models with a narrower halo mass range as mentioned above, with 
higher required cross sections, would also require a larger baryon fraction in DLAs. 

Finally, there may be a possible conflict with observations that limit the luminosity of 
galaxies associated with DLAs. In general, highly luminous galaxies have not been found 
in the proximity of DLAs, although there are cases where they are present: for example, 
a galaxy with a star formation rate of 25M /yr was found coincident with a DLA of very 
high column density by [58]. This star formation rate is just as expected for a 10 12 Mq halo 
dominated by a central galaxy, if as much as half of the baryons in the halo are turning to 
stars in a central galaxy over the age of the universe at z = 2.3. The lack of luminous galaxies 
in most other cases may have various explanations: the star formation rate in a massive halo 
may be distributed over several satellite galaxies, or maybe underestimated observationally 
because of dust absorption. Associated galaxies in massive halos may often be at rather large 
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impact parameters and may sometimes not be recognized for this reason. It is also likely 
that there is a large scatter in the relations of galaxy luminosity and DLA cross section to 
halo mass, and halos with luminous central galaxies may not be the same ones as halos with 
large DLA cross sections. We therefore do not think that our inferred large host halo mass 
for DLAs is necessarily in conflict with the rarity of luminous associated galaxies found so 
far. Precise observational upper limits on the mean luminosity of associated galaxies within 
at least 10 arc seconds of DLAs (corresponding to the expected halo virial radii) would be 
useful. 

6 Conclusions 

The first detection of the cross-correlation of DLAs with the Lya forest absorption on nearby 
lines of sight is presented in this paper. The cross-correlation is well matched by linear theory 
predictions and shows strong evidence for the presence of redshift space distortions. We 
therefore believe that the excess Lya absorption detected around DLAs is explained by the 
average large-scale gravitational mass inflow around the host halos of DLAs that is predicted 
when large-scale structure grows by gravitational evolution. We have used the amplitude of 
the cross-correlation to infer the bias factor of DLAs, related to the mass distribution of the 
host halos. 

The new observational constraint of the bias factor of DLAs, together with the increas- 
ingly accurate determinations of the DLA rate of incidence, lead us to the conclusion that 
a majority of DLAs at z ~ 2 to 3 arise in giant halos of atomic gas clouds around typical 
galaxies in dark matter halos of masses ~ 10 12 Mq. A simple model that works well is that 
DLAs are present in all halos able to accrete photoionized gas (M > 1O 9 M0), with a cross 
section S(M) oc M a , with a ~ 1.1 when ftp = 1, and slightly steeper for larger values of f3p. 
If the mean DLA column density is independent of halo mass, the value a = 1 corresponds 
to the case where a fixed fraction of the halo mass is present in the form of gas in the DLAs, 
so the observations do not require a large increase of the fraction of mass in DLAs with halo 
mass. However, the observations clearly imply that this gas in DLAs must be enormously 
extended in massive halos, with characteristic proper radii ~ 20/ c 1/2 kpc(M/10 12 M & ) a / 2 
which are an increasing fraction of the virial radius as the halo mass increases. 

This leads to a very different picture of high-redshift galaxies compared to their local 
counterparts: typical galaxies at z > 2 are surrounded by atomic gas in the halo with column 
densities comparable to those in present galactic disks, but covering a much larger area and 
containing ~ 10% of all the baryons in the halo, with kinematics that can only be dominated 
by random motions instead of rotation. This is clearly a fundamental fact that must be 
accounted for in any models of galaxy formation. 

The cross-correlation with the Lya forest can be used for other objects, such as quasars 
and metal-line systems, to study their relation with their large-scale environment and infer 
their bias factor. In the case of DLAs, we hope that their bias factor can be measured as a 
function of metallicity as soon as a rough measure of metal abundance can be obtained, which 
is difficult because of the low signal-to-noise of most of the BOSS spectra, the saturation of 
metal lines, and the ionization correction uncertainties. This would set constraints on the 
halo mass - metallicity relation and the correspondence to the mass-metallicity relation for 
galaxies. 
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A Effect of the Mean Transmission Correction on the cross-correlation 

In section 3.3 the Mean Transmission Correction (MTC) has been introduced to reduce the 
variance in the total absorption of each individual spectrum. The MTC consists of renor- 
malizing the transmitted fractions according to equations (3.5) and (3.7). This correction is 
nearly equivalent to requiring the mean perturbation 5pi to be zero when averaged over each 
individual spectrum, except for the fact that the values of F after the MTC are not exactly 
equal to F e , as shown in figure 3. 

The correction was also applied for the study of the Lya correlation function in [13] to 
reduce the variance arising from spectrophotometric errors. As noted there, the correction 
eliminates correlated absorption at all scales, systematically biasing the measured correlation 
function. A method to correct for this systematic bias when fitting a theoretical model to 
the observed correlation function was described in Appendix A of [13], which was tested 
with a set of mocks of Lya spectra. Here, we derive an equivalent correction for the case of 
the cross-correlation of the Lya forest with DLAs, or with any other objects with measured 
redshifts. For the MTC correction used in this paper, we consider only the simple case where 
all pixels have the same weight and instrumental noise is ignored, although the correction 
can be generalized to other cases. We also assume that the average of the corrected Sfi in 
each spectrum is exactly zero. 

Let £d_f(ct, 7r) be the true cross-correlation as a function of the comoving separation 
between the DLA and the Lya pixel in the transverse and parallel directions in redshift 
space, a and ir. We consider a specific Lya forest spectrum that extends from tti to tt2 in 
the parallel separation from the DLA. The corrected cross-correlation ^^(<7, tt), measured 
after the average of 5f has been subtracted from the Lya spectrum, is obtained by simply 
subtracting its value averaged over the range from tt\ to 7T2, 

t) = Zdf(*, tt) — r dn' t DF (a, n') . (A.l) 

The superindex (q) on the corrected cross-correlation indicates that this function depends on 
the specific quasar spectrum and DLA, through the variables tt\ and 112- 

A Lya absorption spectrum measured in a quasar at redshift z q starts and ends at 
absorption redshifts (1 + zi) = Ai(l + z q )/\ a and (I + Z2) = ^2(1 + z q )/X a , with Ai = 1041 A 
and A2 = 1185 A for the analysis presented in this paper. We neglect the fact that quasars 
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with z < 2.4 have a Lya spectrum starting at a higher redshift than Z\ because of the 
wavelength range of the BOSS spectrograph. For a DLA at redshift zd, the values of ir\ 
and 7T2 for this particular spectrum are the comoving separation from zjy to z\ and Z2- The 
correction to apply (second term on the right hand side of equation A.l) to a given pair DLA- 
spectrum is then a function of zd and z q . Defining R as the comoving separation between 
zd and z q , we can express -k\ and TT2 as a function of R and zd- 

In practice we calculate the correction as a function of R only, by considering that 
all our systems are at redshift zd = 2.3, which is the median redshift of the Lya forest 
pixels weighted by their contribution to our cross-correlation measurement, and we neglect 
the changes with redshift except when we measure the cross-correlation at different redshift 
intervals in section 4.8. In order to compute the integral in the right hand side of equation 
A.l we use the linear theory described in section 3. 
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Figure 15. The correction £ cor (a,R) as a function of DLA-quasar radial separation R, for different 
values of a (from top to bottom): (1, 5, 10, 30, 50, 100) comoving hr 1 Mpc. 

The correction to the cross-correlation function is defined as 

t cor (o,R) = Cdf(^tt) - $ F (a,n) = ( % fa! Zdf(?J) • (A.2) 

This correction depends on R, and on the fixed wavelengths Ai and A2, but not on tt. The 
function £ cor is shown in figure 15 as a function of R for several values of a, evaluated at 
z D = 2.3 

The correction is symmetric with respect to a sign change of R — R c , where R c ~ 
265 /i -1 Mpc is the value of R when the DLA redshift is at the center of the forest, i.e., 
7Ti = — TT2- The absolute value of the correction drops sharply for R > 455 fo -1 Mpc or 
R < 75/i _1 Mpc (or equivalently, tt\ > or 712 < 0, respectively), when the DLA redshift is 
outside the Lya forest region of the quasar and the strongest contribution to £ cor near ir = 
is outside the range of the integral in equation (A.2). 

To obtain the average effect on the cross-correlation when measured over our entire 
sample of DLAs and quasar spectra, we define the function p(R\ir) to be the conditional 
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probability of R, given that a Lya forest pixel is present at a comoving parallel separation ir 
from a DLA. This conditional probability depends in a complex way on the redshift distri- 
bution of DLAs and quasars in our sample, but can easily be computed. The final correction 
to the cross-correlation is 

A£(a,7r) = / dRp(R\ir) Scoria, R) . (A.3) 



Conditional probability p(R | it) 
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Figure 16. Conditional probability of R (DLA-quasar radial separation) given 7r (radial separation 
between the Lya pixel and the DLA). p(R\ir) is shown as a function of R for different values of tt, 
computed for our sample of DLAs and quasar spectra. 
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Figure 17. Final correction applied to the theoretical cross-correlation and used when fitting the 
DLA bias. The correction is shown as a function of radial separation 7r for different bins of a (from 
top to bottom): 1 < a < 4 h^ 1 Mpc, 4 < a < 7 h^ 1 Mpc, 7 < a < 1 h^ 1 Mpc, 10 < a < 1 5 h^ 1 Mpc, 
15 < a < 20 hr 1 Mpc, 20 < a < 30 h' 1 Mpc, 30 < a < 40 h' 1 Mpc, and 40 < a < 60 hr 1 Mpc. 

In figure 16 we plot the conditional probability p(R\ir) measured for the distribution 
of DLA and quasar redshifts in our sample, for four different values of n. This conditional 
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probability is approximately invariant if R — tt is fixed, but not exactly so. The quantity 
p(R\tv) is not symmetric under a change of sign of tt, so the correction A£(cr, tt) is also 
asymmetric. This effect is seen in figure 17 which shows the correction we apply as a function 
of tt, for the different bins of a used in our analysis (increasing from top to bottom; see section 
3.4). 

In figure 9 (section 4) we plot the measured cross-correlation function, together with the 
best fit model, with and without the correction. The effect of this correction is clearly not 
crucial for our investigation, but will be more important for studies aiming to measure the 
cross-correlation on larger scales. We expect this to be more important when measuring the 
cross-correlation of quasars and the Lya absorption, because of the larger number of quasars 
and their higher bias factor. 



B Optimal quadratic estimator 

In section 3 our simple method to estimate the cross-correlation of Lya absorption and the 
density of DLAs is described. In this appendix the precise optimal quadratic estimator for 
this cross-correlation is derived, assuming that the fields are Gaussian, and the assumptions 
required to obtain the simple estimator of equation (3.9) are described. 

B.l Definitions and notation 

The estimator of the fluctuation S pi at each pixel i in the whole set of quasar spectra is 

df = 4 -l = S Fi + ^, (B.l) 

Ji Ji 

where is the measured physical flux, /j = C{Fi is the averaged value of fi at redshift Zi, 
and Ni is the noise, which is assumed to have a Gaussian distribution with variance \Nf\. 
The covariance matrix of this estimator is, neglecting any errors from the continuum fit, 

C? 3 = (df df) = fr.(r„) + y%t6(j , (B.2) 

J i 

where £p( r ij) is the correlation function of the Lya fluctuation in the two pixels separated 
by Yij and S K is the Kronecker delta function. This matrix is the same as Cj,- in equation 
3.12, except that the continuum fit errors are not included in the analysis of this appendix. 

The other ingredient that we need to model is the overdensity of detected DLAs. Our 
discussion is completely general and applicable to any set of detected objects that are cross- 
correlated with the Lya forest, and we refer to these detected objects in a specific survey as 
systems. We divide the survey volume in redshift space into iV c small cells of volume V c , such 
that the number hi of systems detected in each cell i is either or 1, and the mean number 
of systems per cell, hi, is much smaller than one. The overdensity in a cell, Shi, is defined 
according to 

hi = hi [1 + S hi ] + ei , (B.3) 
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where ej is the shot-noise, with mean and a variance <r^ = hi. The estimator of the density 
fluctuation from the measured systems in the survey is 

df = £-l = f M +* (B.4) 
hi hi 

and its associated covariance matrix C% is 

C£ = <df df > = fe(r y ) + i-<f , (B.5) 

where £,H( r ij) is the correlation function of detected systems. In general, these detected 
systems may be associated with dark matter halos, and then is the correlation function of 
halos weighted with their probability of yielding a detected system in the survey as a function 
of halo mass. This probability is the mean number of objects per halo when the detected 
systems are galaxies or quasars, or the cross section for the case of DLAs or metal absorption 
lines. 

It is convenient to define a global data vector, 

d = (d F ,d H ), (B.6) 

of dimension Np + N c , where Np is the total number of Lya pixels and N c the number of 
cells. The covariance matrix of the global data vector can be written as 

° = {(C x f C H ) ' (R7) 

where C x is the covariance of the two fields, 

Cf j ={dfdf) = i x {v ij ) , (B.8) 

and £,x( r ij) is the quantity we want to obtain, the cross-correlation function of the Lya 
absorption and the over density of systems. In this appendix the calligraphy letters C are 
used for matrices with dimension (Np + N c ) x (Np + N c ). Generally, we are interested 
in measuring the cross-correlation in some bins in the separation (both the transverse and 
parallel components). Designating the bin number with a subindex a, and the model value 
of the cross-correlation in this bin as £ a , the matrix C x can be expressed as: 

C^ = M^ a , (B.9) 

where Mf- is a matrix with non-zero elements only when the separation between the pixel i 
and the cell j, nj, lies inside the bin r a . 

The derivative of the covariance matrix C with respect to the cross-correlation £ Q in one 
of the bins is 

c -=(a/Jo)' < B10 > 

and the covariance matrix is decomposed as 

C = C t ^ a +M, (B.ll) 

with 

"=(T<M- (B12) 
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B.2 General solution 



An unbiased quadratic estimator can generally be written as ([59], [60]): 



(B.13) 



with 



A a p = Tr[WC a WCp] , 



(B.14) 



and 



Ep = Tr [W C fi W (dd T - M)] 



(B.15) 



where W is any symmetric matrix that we use to weigh our data vector. 

For a Gaussian data vector, the optimal estimator is W = C" 1 , and A is the Fisher 
matrix. Any other weight matrix W still yields an unbiased estimator, but the estimator 
becomes more optimal as W approaches C _1 . 

The expected value of the estimator is 



(L) = (A)~l {E p ) = {A)~l Tr [W C a W(C- AT)] = (A)^ A^ £ y = £ a . (BAG) 



B.3 Approximations 

Inverting covariance matrices is a computationally demanding process. Since we have thou- 
sands of spectra with hundreds of pixels, the matrix C F is in practice impossible to invert, 
not to mention the matrix C. Here, we explain the different approximations that one can do 
in order to get an estimator that is sub-optimal but realistic. 

Noise-dominated regime 

The first approximation is to consider that the contribution of the cross-correlation to the 
covariance is negligible, i.e., we can approximate C" 1 ~ Af^ 1 . The weight matrix is then 




(B.17) 



The product WC a is: 





We can now compute E a and A a p for this weighting: 



E a = Tr [N~ l C, a N~ l (dd T — AT)] = 2Tr [(C F )~ 1 M Q (C H )- 1 d H d^] 



(B.19) 



and 



A a p = Tr [Af^ 1 C t0l M~ 



C, p ] =2Tr[(C F )- 1 M Q (C H )- 1 Mj] 



(B.20) 
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Rare-objects regime 

In this case the covariance matrix of the systems overdensity, C H , will be dominated by the 
shot-noise, and we can do a second approximation for part of the weighting matrix: 

Cf 3 ~ hTH% . (B.21) 

We can now simplify E a even further: 

E a = 2 Tr [(C F )^ 1 M a (C 11 )" 1 d H d£] (B.22) 

N F Np N c N c 

= 2 E E E E( c % lM U cH r k ?dM 

i=0 j=0 k=0 1=0 
N H Np Np 

= 2 EEE^7H F , 

fc=0 i=0 j=0 

where we have used that (C^)" 1 ~ hi df = (hi — hi) equals to one when there is a system 
in the cell and zero otherwise, because the constant hi averages to zero when cross-correlated 
with the Lya field dp. The first summations in k and I are over the N c cells in the volume, 
while the last summation in k is only over those Nh cells with a system on it. 

The matrix A a p can also be simplified in a similar manner: 

A aP = 2 Tr [(C^- 1 M a (C^)" 1 Mj] (B.23) 

N F N F N c N c 

= 2 EEEE( C % 1 ^(^ 1 ^ 

i=0 j=Q k=0 1=0 
N H Np Np 

= 2 EEE^(^^- 

k=0 i=0 j=0 

The last expressions can be rewritten in a more compact form by defining a new matrix 
M' a , which has dimension Njj x Np (instead of N c x Np for the matrix M a ) because it is 
defined only in the cells where a system is detected, 

E a = 2 Tr [(C F )~ 1 M' a 1 1 d%] , (B.24) 

and 

A a/s = 2TT[(C F )~ 1 M' a lM' p T ] . (B.25) 

We have included the identity matrix 1 and vector 1, both with dimension Nh, to ease the 
reference to the previous expressions. 

Independent-spectra regime 

Pixels of two different spectra are only weakly correlated, especially in the limit of low quasar 
density where the typical separation between quasar lines of sight is large. In this limit, an 
additional approximation can be made for the weighting matrix in which only pixels from 
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the same spectrum are correlated. Then, the correlation matrix C F and its inverse are block 
diagonal, with a block C F for each quasar: 



/Cf .. 


. . 


• \ 


C F .. 


. . 


. 


0.. 


'.cf'. 


. 


\ o o .. 


. . 





(B.26) 



Instead of having to invert a matrix of size Np x Np (with Np ~ 5 x 10 7 ) we need only 
N q matrices of size N p x N p , where N q ~ 10 5 is the number of spectra and N p ~ 500 is the 
number of pixels in a typical spectrum. 



Independent-pixels regime 

Our final approximation is that each pixel in the spectrum is independent. In this case the 
matrix C F is no longer block diagonal, but diagonal: 



C ij = &j + ~yr~ °ij ~ a Ti °ij > (B.27) 
/ i 



where a\ i = a Fi + (N) / ff is the total variance in the pixel, including the Lya intrinsic 
variance. 

The expressions for E a and A a p are now equal to 

N H Np N F N H N F 

E « = 2 E E E MjkVW d f = 2 E E M % *>* d f > (b.28) 

fc=0 i=0 j=0 fc=0 «=0 

and 

JV H jV> N F N H Np 

= 2 E E E ^(c^ 1 M l = 2 E E ^ M £ > (b.29) 

k=0 i=0 j=0 fc=0 i=0 

where we have defined the pixel weight Wi = . 

We now see that the estimator in equation 3.9 corresponds to this case if we use only 
the diagonal part of A a p\ 

N H N F N H 

E E M tk w i 4 EE^rff 

t A — 1 1? fe=0i=0 fc=0«€a / c on-i 

= A aa E a = = — . (B.30) 

E E M ik w i M ik EE«i 

fc=0i=0 fc=0i£ct 



If wide redshift bins are used with a substantial evolution of the signal being present 
within the bin, the estimator needs to be modified. In particular, equation B.9 should now 

X i3 =M 3f „(l±£)\ (B.31) 
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where we have assumed that the evolution of the cross-correlation is described by an am- 
plitude varying as a power-law with the scale factor with an index 7. The simple estimator 
used in Section 3 needs to be modified as well, by multiplying the weights in equation 3.10 
by the same factor. 

We have not applied this correction for the weights, because the cross-correlation signal 
does not evolve strongly with redshift. In an Einstein-de Sitter universe, or at sufficiently 
high redshift to make the dark energy component negligible, the growth factor is proportional 
to the scale factor, and the value of 7 is simply related to the evolution of the bias parameter: 



I + ZijY 1 _ ( b(zij)\ fl + Zij 



1 + z a J V H z a) J V 1 + z - 



(B.32) 



From the evolution of the line-of-sight Lya forest power spectrum measured in [36], we can 
assume that the amplitude of the Lya clustering evolves with 7^ = 3.8, implying that bp 
evolves as (1 + z) 2 ' 9 . Assuming that the DLA bias is constant over the redshift range, we 
obtain a value of jx = 0.9 for the cross-correlation. This leads to a very small change over 
the narrow redshift range of our DLA sample. 
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